#define dataselector_cxx
// The class definition in dataselector.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("dataselector.C")
// Root > T->Process("dataselector.C","some options")
// Root > T->Process("dataselector.C+")
//

#include "dataselector.h"
#include <TH2.h>
#include <TMath.h>
#include <TStyle.h>
#include "PlotQCD.C"

void dataselector::Begin(TTree * /*tree*/)
{
   // The Begin() function is called at the start of the query.
   // When running with PROOF Begin() is only called on the client.
   // The tree argument is deprecated (on PROOF 0 is passed).

   cout << "Begin" << endl;

   TString option = GetOption();

   // setting runtime options
   ParseOptions(option);

   // fake efficiency summary table 
   d0sig_ctl      = 0;
   num_jets       = -1;
   num_mu_tight   = 0.;
   num_mu_loose   = 0.;
   eff_fake       = 0.;
   eff_fake_error = 0.;

   if (is_data || is_mc_mix || is_mc_zjets || is_mc_wjets || is_mc_ttbar ) {
      tree_eff_fake    = new TTree(sample + "_eff_fake","");
      b_d0sig_ctl      = tree_eff_fake->Branch("d0sig_ctl", &d0sig_ctl, "d0sig_ctl/D");
      b_num_jets       = tree_eff_fake->Branch("num_jets", &num_jets, "num_jets/I");
      b_num_mu_tight   = tree_eff_fake->Branch("num_mu_tight", &num_mu_tight, "num_mu_tight/D");
      b_num_mu_loose   = tree_eff_fake->Branch("num_mu_loose", &num_mu_loose, "num_mu_loose/D");
      b_eff_fake       = tree_eff_fake->Branch("eff_fake", &eff_fake, "eff_fake/D");
      b_eff_fake_error = tree_eff_fake->Branch("eff_fake_error", &eff_fake_error, "eff_fake_error/D");
   }

   // fake efficiency histograms 
   for (int j=0; j<hist_array_size; j++) {

      TString ht_ext = getHistNameExt(j);

      hEff[j] = new TH1D(sample + "_eff_l2t" + ht_ext ,"", 12, 0, 24.);
   }
}

void dataselector::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   cout << "SlaveBegin" << endl;

   TString option = GetOption();

   // setting runtime options
   ParseOptions(option);

   // initialize the output histograms
   N_event = 0;

   // histogram initialization
   InitHistograms();

   // pack output from slave into outputlist
   PackOutputList();
}

Bool_t dataselector::Process(Long64_t entry)
{
   // The Process() function is called for each entry in the tree (or possibly
   // keyed object in the case of PROOF) to be processed. The entry argument
   // specifies which entry in the currently loaded tree is to be processed.
   // It can be passed to either dataselector::GetEntry() or TBranch::GetEntry()
   // to read either all or the required parts of the data. When processing
   // keyed objects with PROOF, the object is already loaded and is available
   // via the fObject pointer.
   //
   // This function should contain the "body" of the analysis. It can contain
   // simple or elaborate selection criteria, run algorithms on the data
   // of the event and typically fill histograms.
   //
   // The processing can be stopped by calling Abort().
   //
   // Use fStatus to set the return value of TTree::Process().
   //
   // The return value is currently not used.

   // Get the entire event
   // GetEntry(entry);

   if ( N_event > 0 && N_event % 1000000 == 0 ) cout << "processed " << N_event << " events" << endl;

   N_event += 1;

   /* == NOTE == */
   /* applying event-level cuts based on the TopObservationSelection */
   /* == NOTE == */

   // 1. good run list on data
   b_grlevent->GetEntry(entry);
   if ( is_data && (! grlevent) ) return kFALSE;

   // 2. with one good muon
   //if ( ! GoodLeptonCut ) return kFALSE;
   //b_GoodLeptonCut->GetEntry(entry);
   b_GoodMuonCut->GetEntry(entry);
   //if ( ! (GoodMuonCut && GoodLeptonCut) ) return kFALSE;

   // with no more than 1 good lepton
   b_GoodMu_N->GetEntry(entry);
   b_GoodEl_N->GetEntry(entry);
   if ( ! (GoodMu_N + GoodEl_N < 2) ) return kFALSE;

   // 3. HFORCut (??)
   b_HFORCut->GetEntry(entry);
   if ( ! HFORCut ) return kFALSE;

   // 4. trigger selection and trigger-match dep. on runs
   b_TriggerCut->GetEntry(entry);
   if ( ! TriggerCut ) return kFALSE;
   //b_TriggerMatchCut->GetEntry(entry);
   //if ( ! (TriggerCut && TriggerMatchCut) ) return kFALSE;

   // 5. missing Et > 20 GeV 
   b_EtMissCut->GetEntry(entry);
   if ( ! EtMissCut ) return kFALSE;

   // 6. having primiry vertex 
   b_NonColBkgCut->GetEntry(entry);
   if ( ! NonColBkgCut ) return kFALSE;

   // 7. no bad jets 
   b_CleaningCut->GetEntry(entry);
   if ( ! CleaningCut ) return kFALSE;

   // 8. lepton overlap
   b_emuOverlapCut->GetEntry(entry);
   if ( ! emuOverlapCut ) return kFALSE;

   // 9. number of Jets with energy > 25 GeV 
   b_GoodJet25_N->GetEntry(entry);
   if ( GoodJet25_N < 1 ) return kFALSE;

   // 10. btagging cuts (optional)
   if ( btagging_cut ) {
       b_BtagCut->GetEntry(entry);
       if ( ! BtagCut ) return kFALSE;
   }

   // 11. triangule cuts (optional) - triangular cut treated in particle level on loose Muons
   b_TriangularCut->GetEntry(entry);

   /* == NOTE == */
   /* distinguishing loose and tight muons in the selected events */
   /* == NOTE == */
   b_GoodMu_dR_jet->GetEntry(entry);
   b_GoodMu_etcone30->GetEntry(entry);
   b_GoodMu_ptcone30->GetEntry(entry);
   b_GoodMu_d0_sig_exPV->GetEntry(entry);
   b_GoodMu_pt->GetEntry(entry);
   b_GoodMu_phi->GetEntry(entry);

   b_topMET_et->GetEntry(entry);
   b_topMET_etx->GetEntry(entry);
   b_topMET_ety->GetEntry(entry);

   b_MWtrans->GetEntry(entry);
   b_topMass_GoodJet->GetEntry(entry);
   b_NBjets25->GetEntry(entry);

   int jet_idx = GoodJet25_N;
   if ( jet_idx > max_jet_multi ) jet_idx = max_jet_multi;

   Double_t my_weight = 1.0;

   // presuming the MC NTuple is normalized to 100 pb^-1
   if ( is_mc ) {
       b_mcevt_weight->GetEntry(entry);
       b_weight->GetEntry(entry);
       my_weight = (mcevt_weight * weight) / 100.;
   }

   if ( is_data || is_mc_qcd || is_mc_mix || is_mc_ttbar || is_mc_wjets || is_mc_zjets ) {

      b_LooseMu_N->GetEntry(entry);
      b_LooseMu_dR_jet->GetEntry(entry);
      b_LooseMu_etcone30->GetEntry(entry);
      b_LooseMu_ptcone30->GetEntry(entry);
      b_LooseMu_d0_sig_exPV->GetEntry(entry);
      b_LooseMu_pt->GetEntry(entry);
      b_LooseMu_phi->GetEntry(entry);

      for(int i=0; i<LooseMu_N; i++) {
 
         // pass loose selection
         if ( LooseMu_dR_jet->at(i) > 0.4 ) {
 
             // calculate W transverse mass with loose muons 
             Double_t my_mwt = sqrt( pow( topMET_et + LooseMu_pt->at(i), 2) - pow( topMET_etx + TMath::Cos(LooseMu_phi->at(i))*LooseMu_pt->at(i), 2) - pow( topMET_ety + TMath::Sin(LooseMu_phi->at(i))*LooseMu_pt->at(i), 2) );
 
             // triangle cut on muon level
             if ( triangle_cut && ( my_mwt + topMET_et < 60000. ) ) continue;
 
             hLoose[0]->Fill(LooseMu_d0_sig_exPV->at(i), my_weight);
             hLoose[jet_idx]->Fill(LooseMu_d0_sig_exPV->at(i), my_weight);
  
             if ( LooseMu_d0_sig_exPV->at(i) > d0_sig_t ) {
                 hNLoose[0]->Fill(1, my_weight);
                 hNLoose[jet_idx]->Fill(1, my_weight);
             }
 
             hMWt_loose[0]->Fill(my_mwt, my_weight);
             hMWt_loose[jet_idx]->Fill(my_mwt, my_weight);
             hNJets_loose[0]->Fill(GoodJet25_N, my_weight);
             hNJets_loose[jet_idx]->Fill(GoodJet25_N, my_weight);
             hNBjets25_loose[0]->Fill(NBjets25, my_weight);
             hNBjets25_loose[jet_idx]->Fill(NBjets25, my_weight);
 
             if (topMET_et > 0) {
                 hMET_loose[0]->Fill(topMET_et, my_weight);
                 hMET_loose[jet_idx]->Fill(topMET_et, my_weight);
             }
 
             if (topMass_GoodJet > 0) { 
                 hMTop_loose[0]->Fill(topMass_GoodJet, my_weight);
                 hMTop_loose[jet_idx]->Fill(topMass_GoodJet, my_weight);
             }
         }
      }

      if ( ! GoodMuonCut ) return kFALSE;
  
      for(int i=0; i<GoodMu_N; i++) {
 
         // pass tight selection
         if ( GoodMu_dR_jet->at(i) > 0.4 && GoodMu_etcone30->at(i) < 4000. && GoodMu_ptcone30->at(i) < 4000. ) {
 
             // calculate W transverse mass with tight muons 
             Double_t my_mwt = sqrt( pow( topMET_et + GoodMu_pt->at(i), 2) - pow( topMET_etx + TMath::Cos(GoodMu_phi->at(i))*GoodMu_pt->at(i), 2) - pow( topMET_ety + TMath::Sin(GoodMu_phi->at(i))*GoodMu_pt->at(i), 2) );
 
             // triangle cut on muon level
             if ( triangle_cut && ( my_mwt + topMET_et < 60000. ) ) continue;
 
             hTight[0]->Fill(GoodMu_d0_sig_exPV->at(i), my_weight);
             hTight[jet_idx]->Fill(GoodMu_d0_sig_exPV->at(i), my_weight);
 
             if ( GoodMu_d0_sig_exPV->at(i) > d0_sig_t ) {
                 hNTight[0]->Fill(1, my_weight);
                 hNTight[jet_idx]->Fill(1, my_weight);
             }
 
             // calculate observables with tight selection
             hMWt_tight[0]->Fill(my_mwt, my_weight);
             hMWt_tight[jet_idx]->Fill(my_mwt, my_weight);
 
             hNJets_tight[0]->Fill(GoodJet25_N, my_weight);
             hNJets_tight[jet_idx]->Fill(GoodJet25_N, my_weight);
             hNBjets25_tight[0]->Fill(NBjets25, my_weight);
             hNBjets25_tight[jet_idx]->Fill(NBjets25, my_weight);
 
             if (topMET_et > 0) {
                 hMET_tight[0]->Fill(topMET_et, my_weight);
                 hMET_tight[jet_idx]->Fill(topMET_et, my_weight);
             }
 
             if (topMass_GoodJet > 0) { 
                 hMTop_tight[0]->Fill(topMass_GoodJet, my_weight);
                 hMTop_tight[jet_idx]->Fill(topMass_GoodJet, my_weight);
             }
         }
      }
   }

   // observables available in NTuple
   //  - presuming the event selection for top analysis has been done in NTuple making
   //  - applying only the missing cut in the selection implemented in TopInputs
   //  - the histograms filled below should be identical to the histograms filled with tight selection on loose muons 
   if ( GoodMuonCut && GoodMu_N > 0 && GoodMu_dR_jet->at(0) > 0.4 && GoodMu_etcone30->at(0) < 4000 && GoodMu_ptcone30->at(0) < 4000 ) {

       // triangle cut
       if ( triangle_cut && (! TriangularCut) ) return kFALSE;

       if ( MWtrans > 0 ) {
           hMWt[0]->Fill(MWtrans, my_weight);
           hMWt[jet_idx]->Fill(MWtrans, my_weight);
       }

       if ( topMass_GoodJet > 0 ) { 
           hMTop[0]->Fill(topMass_GoodJet, my_weight);
           hMTop[jet_idx]->Fill(topMass_GoodJet, my_weight);
       }

       if ( topMET_et > 0 ) {
           hMET[0]->Fill(topMET_et, my_weight);
           hMET[jet_idx]->Fill(topMET_et, my_weight);
       }
 
       hNJets[0]->Fill(GoodJet25_N, my_weight);
       hNJets[jet_idx]->Fill(GoodJet25_N, my_weight);

       hNBjets25[0]->Fill(NBjets25, my_weight);
       hNBjets25[jet_idx]->Fill(NBjets25, my_weight);
   }

   return kTRUE;
}

void dataselector::SlaveTerminate()
{
   // The SlaveTerminate() function is called after all entries or objects
   // have been processed. When running with PROOF SlaveTerminate() is called
   // on each slave server.

   cout << "SlaveTerminate" << endl;
}

void dataselector::Terminate()
{
   // The Terminate() function is the last function to be called during
   // a query. It always runs on the client, it can be used to present
   // the results graphically or save the results to file.

   cout << "Terminate" << endl;

   // unpack output from slave into outputlist
   UnpackOutputList();

   if ( is_data || is_mc_qcd || is_mc_mix || is_mc_ttbar || is_mc_wjets || is_mc_zjets ) {
       
      // filling up the loose2tight efficiency histogram
      for (int j=0; j < hist_array_size; j++) {
 
         for (Int_t i=0; i < hLoose[j]->GetNbinsX() + 1; i++) {
         
              Double_t iloose = hLoose[j]->GetBinContent(i);
              Double_t itight = hTight[j]->GetBinContent(i);
 
              cout << "bin " << i << ":" << itight << " " << iloose << endl;
         
              Double_t eff = ( iloose <= 0.) ? 0.:itight/iloose;
         
              Double_t dN = iloose - itight;
         
              Double_t stat_err = 0.0;
         
              if ( eff != 0. ) {
                  Double_t stat_err = TMath::Sqrt( 
                                        TMath::Power( (iloose + itight) / TMath::Power(iloose,2) , 2) * itight +
                                        TMath::Power( itight / TMath::Power(iloose,2), 2) * dN
                                      );
              }
         
              hEff[j]->SetBinContent(i,eff);
              hEff[j]->SetBinError(i,stat_err);
         }
      }
   }

   // count the number of loose and tight muons and calculate the fake efficiency in various samples
   if (is_data || is_mc_mix || is_mc_zjets || is_mc_wjets || is_mc_ttbar ) {

       for (int j=0; j < hist_array_size; j++) {

           //Double_t nloose = hNLoose[j]->GetEntries();
           //Double_t ntight = hNTight[j]->GetEntries();

           // scale to data luminosity before counting
           hNLoose[j]->Scale(h_scale);
           hNTight[j]->Scale(h_scale);

           Double_t nloose = hNLoose[j]->Integral();
           Double_t ntight = hNTight[j]->Integral();

           cout << "j=" << j << ":" << h_scale << " " << nloose << " " << ntight << endl;

           // calculate the fake efficiency
           eff_fake = ( nloose <= 0) ? 0.:( ntight / nloose );
         
           // calculate the statistic error of the fake efficiency
           Double_t dN = nloose - ntight;
         
           eff_fake_error = 0.;
         
           if ( eff_fake != 0. ) {
                eff_fake_error = TMath::Sqrt(
                                    TMath::Power( (nloose + ntight) / TMath::Power(nloose,2) , 2) * ntight +
                                    TMath::Power( ntight / TMath::Power(nloose,2), 2) * dN
                                 );
           }

           cout << "jet multiplicity:" << j              << endl;
           cout << "QCD MM - Ef     :" << eff_fake       << endl;
           cout << "QCD MM - Ef_err :" << eff_fake_error << endl;

           num_jets     = j;
           d0sig_ctl    = d0_sig_t;
           num_mu_tight = ntight;
           num_mu_loose = nloose;
           tree_eff_fake->Fill();

           // call data-driven QCD MM only on data sample
           if ( is_data ) {

               TString ht_ext = "_j";
               ht_ext += j;
            
               hMWt_qcd_estat[j] = PlotQCD(hMWt_loose[j], hMWt_tight[j], eff_fake, eff_fake_error, 4, 0, 0);
               hMWt_qcd_estat[j]->SetNameTitle(sample+"_mwt_qcd_estat"+ht_ext,"");
               hMWt_qcd_esyst[j] = PlotQCD(hMWt_loose[j], hMWt_tight[j], eff_fake, eff_fake_error, 4, 0, 1);
               hMWt_qcd_esyst[j]->SetNameTitle(sample+"_mwt_qcd_esyst"+ht_ext,"");
               hMWt_qcd[j] = PlotQCD(hMWt_loose[j], hMWt_tight[j], eff_fake, eff_fake_error, 4, 0, 2);
               hMWt_qcd[j]->SetNameTitle(sample+"_mwt_qcd"+ht_ext,"");
         
               hMTop_qcd_estat[j] = PlotQCD(hMTop_loose[j], hMTop_tight[j], eff_fake, eff_fake_error, 4, 0, 0);
               hMTop_qcd_estat[j]->SetNameTitle(sample+"_mtop_qcd_estat"+ht_ext,"");
               hMTop_qcd_esyst[j] = PlotQCD(hMTop_loose[j], hMTop_tight[j], eff_fake, eff_fake_error, 4, 0, 1);
               hMTop_qcd_esyst[j]->SetNameTitle(sample+"_mtop_qcd_esyst"+ht_ext,"");
               hMTop_qcd[j] = PlotQCD(hMTop_loose[j], hMTop_tight[j], eff_fake, eff_fake_error, 4, 0, 2);
               hMTop_qcd[j]->SetNameTitle(sample+"_mtop_qcd"+ht_ext,"");
         
               hMET_qcd_estat[j] = PlotQCD(hMET_loose[j], hMET_tight[j], eff_fake, eff_fake_error, 4, 0, 0);
               hMET_qcd_estat[j]->SetNameTitle(sample+"_met_qcd_estat"+ht_ext,"");
               hMET_qcd_esyst[j] = PlotQCD(hMET_loose[j], hMET_tight[j], eff_fake, eff_fake_error, 4, 0, 1);
               hMET_qcd_esyst[j]->SetNameTitle(sample+"_met_qcd_esyst"+ht_ext,"");
               hMET_qcd[j] = PlotQCD(hMET_loose[j], hMET_tight[j], eff_fake, eff_fake_error, 4, 0, 2);
               hMET_qcd[j]->SetNameTitle(sample+"_met_qcd"+ht_ext,"");
         
               hNJets_qcd_estat[j] = PlotQCD(hNJets_loose[j], hNJets_tight[j], eff_fake, eff_fake_error, 4, 0, 0);
               hNJets_qcd_estat[j]->SetNameTitle(sample+"_njets_qcd_estat"+ht_ext,"");
               hNJets_qcd_esyst[j] = PlotQCD(hNJets_loose[j], hNJets_tight[j], eff_fake, eff_fake_error, 4, 0, 1);
               hNJets_qcd_esyst[j]->SetNameTitle(sample+"_njets_qcd_esyst"+ht_ext,"");
               hNJets_qcd[j] = PlotQCD(hNJets_loose[j], hNJets_tight[j], eff_fake, eff_fake_error, 4, 0, 2);
               hNJets_qcd[j]->SetNameTitle(sample+"_njets_qcd"+ht_ext,"");
            
               hNBjets25_qcd_estat[j] = PlotQCD(hNBjets25_loose[j], hNBjets25_tight[j], eff_fake, eff_fake_error, 4, 0, 0);
               hNBjets25_qcd_estat[j]->SetNameTitle(sample+"_nbjets25_qcd_estat"+ht_ext,"");
               hNBjets25_qcd_esyst[j] = PlotQCD(hNBjets25_loose[j], hNBjets25_tight[j], eff_fake, eff_fake_error, 4, 0, 1);
               hNBjets25_qcd_esyst[j]->SetNameTitle(sample+"_nbjets25_qcd_esyst"+ht_ext,"");
               hNBjets25_qcd[j] = PlotQCD(hNBjets25_loose[j], hNBjets25_tight[j], eff_fake, eff_fake_error, 4, 0, 2);
               hNBjets25_qcd[j]->SetNameTitle(sample+"_nbjets25_qcd"+ht_ext,"");
           }
       }
   }

   // make up the output file name 
   ostringstream stm;
   stm << "qcd_mm_d0sig_" << d0_sig_t;
   if ( btagging_cut ) stm << "_btag";
   if ( triangle_cut ) stm << "_tria";
   stm << ".root";
   output_fname = TString( stm.str() );

   // write out histograms 
   //TFile *f = new TFile(sample + ".root","RECREATE");
   TFile *f = new TFile(output_fname,"UPDATE");

   for (int j=0; j<hist_array_size; j++) {

      // observable histograms
      hMWt[j]->Scale(h_scale);      hMWt[j]->Write();
      hMTop[j]->Scale(h_scale);     hMTop[j]->Write();
      hMET[j]->Scale(h_scale);      hMET[j]->Write();
      hNJets[j]->Scale(h_scale);    hNJets[j]->Write();
      hNBjets25[j]->Scale(h_scale); hNBjets25[j]->Write();

      if ( is_data || is_mc_qcd || is_mc_mix || is_mc_ttbar || is_mc_wjets || is_mc_zjets ) {

          // histograms for QCD MM studies
          hNLoose[j]->Write();
          hNTight[j]->Write();
          hTight[j]->Write();
          hLoose[j]->Write();
          hTight[j]->Write();
          hEff[j]->Write();

          hMWt_tight[j]->Scale(h_scale); hMWt_tight[j]->Write();
          hMWt_loose[j]->Scale(h_scale); hMWt_loose[j]->Write();
  
          hMTop_tight[j]->Scale(h_scale); hMTop_tight[j]->Write();
          hMTop_loose[j]->Scale(h_scale); hMTop_loose[j]->Write();
  
          hMET_tight[j]->Scale(h_scale); hMET_tight[j]->Write();
          hMET_loose[j]->Scale(h_scale); hMET_loose[j]->Write();
  
          hNJets_tight[j]->Scale(h_scale); hNJets_tight[j]->Write();
          hNJets_loose[j]->Scale(h_scale); hNJets_loose[j]->Write();
  
          hNBjets25_tight[j]->Scale(h_scale); hNBjets25_tight[j]->Write();
          hNBjets25_loose[j]->Scale(h_scale); hNBjets25_loose[j]->Write();
      }

      if ( is_data ) {
          // data driven QCD estimation histograms
          hMWt_qcd[j]->Scale(h_scale);   hMWt_qcd[j]->Write();
          hMWt_qcd_estat[j]->Scale(h_scale);   hMWt_qcd_estat[j]->Write();
          hMWt_qcd_esyst[j]->Scale(h_scale);   hMWt_qcd_esyst[j]->Write();

          hMTop_qcd[j]->Scale(h_scale);  hMTop_qcd[j]->Write();
          hMTop_qcd_estat[j]->Scale(h_scale);  hMTop_qcd_estat[j]->Write();
          hMTop_qcd_esyst[j]->Scale(h_scale);  hMTop_qcd_esyst[j]->Write();

          hMET_qcd[j]->Scale(h_scale);   hMET_qcd[j]->Write();
          hMET_qcd_estat[j]->Scale(h_scale);   hMET_qcd_estat[j]->Write();
          hMET_qcd_esyst[j]->Scale(h_scale);   hMET_qcd_esyst[j]->Write();

          hNJets_qcd[j]->Scale(h_scale); hNJets_qcd[j]->Write();
          hNJets_qcd_estat[j]->Scale(h_scale); hNJets_qcd_estat[j]->Write();
          hNJets_qcd_esyst[j]->Scale(h_scale); hNJets_qcd_esyst[j]->Write();

          hNBjets25_qcd[j]->Scale(h_scale); hNBjets25_qcd[j]->Write();
          hNBjets25_qcd_estat[j]->Scale(h_scale); hNBjets25_qcd_estat[j]->Write();
          hNBjets25_qcd_esyst[j]->Scale(h_scale); hNBjets25_qcd_esyst[j]->Write();
      }
   }
   
   if (is_data || is_mc_mix || is_mc_zjets || is_mc_wjets || is_mc_ttbar ) tree_eff_fake->Write();

   // write out a skimmed tree selected events 
   //fChainSelect->Write();

   f->Close();

   // clean all objects in the output list
   CleanOutputList();

   for (int j=0; j<hist_array_size; j++) {
      delete hLoose[j];
      delete hTight[j];
      delete hEff[j];

      delete hMWt[j];
      delete hMWt_tight[j];
      delete hMWt_loose[j];

      delete hMTop[j];
      delete hMTop_tight[j];
      delete hMTop_loose[j];

      delete hMET[j];
      delete hMET_tight[j];
      delete hMET_loose[j];

      delete hNJets[j];
      delete hNJets_tight[j];
      delete hNJets_loose[j];

      delete hNBjets25[j];
      delete hNBjets25_tight[j];
      delete hNBjets25_loose[j];

      if ( is_data ) {
          delete hMWt_qcd[j];      delete hMWt_qcd_estat[j];      delete hMWt_qcd_esyst[j];
          delete hMTop_qcd[j];     delete hMTop_qcd_estat[j];     delete hMTop_qcd_esyst[j];
          delete hMET_qcd[j];      delete hMET_qcd_estat[j];      delete hMET_qcd_esyst[j];
          delete hNJets_qcd[j];    delete hNJets_qcd_estat[j];    delete hNJets_qcd_esyst[j];
          delete hNBjets25_qcd[j]; delete hNBjets25_qcd_estat[j]; delete hNBjets25_qcd_esyst[j];
      }
   }

   if (is_data || is_mc_mix || is_mc_zjets || is_mc_wjets || is_mc_ttbar ) delete tree_eff_fake;
}
